Optimal  LES  and  a  New  Near- Wall  Model  for  Application  to 
High-Reynolds-Number  Air  Foils 

Final  Report  for  Contract  F9550-07- 1-041 1 

R.  D.  Moser 

Institute  for  Computational  Engineering  and  Sciences  and 
Department  of  Mechanical  Engineering 
University  of  Texas  at  Austin 
Austin,  TX  78712 


1  Motivation  and  Objectives 


One  of  the  most  promising  techniques  for  the  prediction  of  turbulent  flows  is  Large  Eddy  Simu¬ 
lation  (LES),  in  which  the  largest  scales  of  turbulent  fluid  motion  are  simulated  while  the  effect 
of  the  smaller  unresolved  scales  are  modeled.  However,  one  of  the  major  obstacles  to  the  use  of 
LES  in  technologically  important  turbulent  flows,  such  as  the  flow  over  an  airfoil,  is  the  modeling 
of  the  near-wall  turbulence.  Current  LES  modeling  approaches  require  that  either  the  near-wall 
turbulence  be  adequately  resolved  (at  unacceptable  expense  for  large  Reynolds  numbers),  or  that 
an  LES  wall-model  be  used,  which  to  date  has  not  provided  accurate  results  in  relatively  complex 
flows  (e.g.  an  airfoil  near  stall). 

A  new  modeling  approach  for  wall-bounded  turbulence  has  been  developed,  which  promises  to 
address  the  wall-modeling  problem  described  above.  It  is  based  on  a  formal  filtering  of  the  tur¬ 
bulence  and  the  boundary,  a  determination  of  the  wall  stress  to  minimize  the  “leakage”  of  kinetic 
energy  and  momentum  from  the  flow  domain,  and  the  use  of  the  optimal  LES  formulation.  These 
three  ingredients  have  been  shown  to  produce  accurate  LES  results  in  a  turbulent  channel.  Further, 
each  of  these  modeling  formulations  is  generally  applicable,  so  they  can  be  expected  to  provide 
a  basis  for  successful  wall  modeling  in  more  complex  flow  situations,  such  as  the  flow  over  an 
airfoil. 

The  objective  of  this  research  was  thus  to  further  refine  and  validate  this  wall-modeling  approach, 
and  ultimately  apply  it  to  the  flow  over  an  airfoil.  To  accomplish  this,  we  pursued  a  number  of 
intermediate  objectives: 


•  Completion  and  testing  of  the  optimal  LES  (OLES)  jfinite  volume  formulation  using  only 
theoretical  inputs.  Such  theory-based  OLES  is  necessary  so  that  the  models  will  be  predictive 
(i.e.  will  not  need  empirical  inputs). 

•  Refinement  and  generalization  of  an  approximation  of  the  multi-point  correlations  in  wall- 
bounded  turbulence.  A  formulation  for  the  anisotropy  and  inhomogeneity  of  the  correlations 
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is  needed  to  apply  the  theoretical  formulation  described  above  for  strongly  inhomogeneous 
turbulence  as  occurs  nears  walls. 

•  Generalization  of  the  wall-bounded  formulation  used  in  [6]  to  the  finite  volume  OLES  for¬ 
mulation,  which  is  appropriate  for  use  in  general  geometries. 

•  Application  of  the  above  OLES  generalizations  to  general  wall-bounded  flows. 


2  Background  and  Approach 


It  is  well  known  that  our  ability  to  predict  complex  fluid  flows  is  limited  by  the  need  to  model 
the  effects  of  turbulence.  For  example,  Reynolds  averaged  (RANS)  turbulence  models  have  dif¬ 
ficulties  in  flows  exhibiting  large  separation  and  vortex  shedding.  Thus,  in  the  technologically 
critical  problem  of  flow  over  an  airfoil,  there  remains  a  great  need  to  improve  turbulence  predictive 
techniques,  especially  at  large  angle  of  attack.  Large  eddy  simulation  (LES)  is  increasingly  being 
employed  to  avoid  the  weaknesses  of  RANS  models,  and  to  provide  more  information  than  RANS 
computations  are  able  to  provide  (e.g.  for  aeroacoustic  phenomena  or  unsteady  forces).  However, 
LES  remains  prohibitively  expensive  for  high  Reynolds  number  applications,  including  airfoils, 
largely  due  to  the  treatment  of  walls  [31,  23]. 

Recently,  a  European  consortium  undertook  to  evaluate  current  capabilities  of  LES  for  application 
to  airfoils  [23].  A  high  lift  airfoil  at  angle  of  attack  sufficient  to  produce  a  mild  separation  on 
the  suction  side  toward  the  trailing  edge  was  simulated  by  a  variety  of  groups,  using  an  array  of 
different  LES  modeling  techniques.  In  this  study,  it  was  found  that  with  current  LES  techniques,  it 
was  necessary  to  resolve  the  near- wall  turbulence  to  obtain  accurate  results.  LES  performed  with 
wall  function  representations  of  the  near-wall  region  were  not  satisfactory. 

The  cost  implications  of  resolving  the  near  wall  layer  are  severe,  with  classical  estimates  for  sim¬ 
ulating  the  boundary  layer  (as  on  an  airfoil)  in  LES  yielding  costs  that  scale  with  Re°cb  if  the 
near-wall  layer  is  not  resolved,  versus  Re2A  if  the  near-wall  layer  must  be  resolved,  as  suggested 
above  [34,  10].  The  reason  for  the  difference  is  that  resolution  of  the  near- wall  layer  requires  grid 
sizes  that  scale  with  the  viscous  wall  unit  {v/uT,  where  uT  is  the  friction  velocity,  defined  in  terms 
of  the  mean  wall  shear  stress  u 2  =  tw/ p),  which  get  small  relative  to  the  boundary  layer  thick¬ 
ness  approximately  like  Re~l.  As  a  result,  such  wall-resolving  LES  have  only  been  performed  for 
airfoils  at  moderate  Reynolds  numbers,  with  very  narrow  spanwise  domains.  For  example,  [22] 
used  a  spanwise  domain  size  of  just  1.2%  of  chord  to  simulate  the  A- Airfoil  from  Aerospatiale  at 
a  chord  Reynolds  number  of  2.1  x  106.  Wall-resolved  LES  of  an  actual  three  dimensional  wing  is 
currently  out  of  the  question. 

It  is  clear  that  an  LES  wall  modeling  approach  is  needed  that  does  not  require  the  resolution  of 
the  near-wall  layer,  and  a  number  of  research  efforts  have  been  directed  at  this  problem,  with  only 
limited  success.  For  an  account  of  these  efforts,  see  the  recent  review  by  [31],  as  well  as  several 
more  recent  papers  [41,  37]  The  approaches  generally  pose  stress  boundary  conditions  for  the  LES 
equations,  and  a  number  of  modeling  approaches  to  determine  the  required  boundary  stresses  have 
been  developed.  These  include  the  correlation  of  the  wall  stresses  with  large-scale  velocities  in  the 
interior  [32],  a  boundary  layer  representation  of  the  wall  layer  [41]  and  the  use  of  optimal  control 
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to  ensure  that  known  statistical  properties  of  the  turbulence  are  reproduced  [37,  27]. 

Under  a  previous  AFOSR/NSF  jointly  funded  project,  a  new  approach  to  the  wall-modeling  prob¬ 
lem  was  shown  to  yield  remarkably  good  results  in  turbulent  channel  flow  [11,  25].  This  approach 
is  based  on  the  observation  that  in  an  LES,  it  is  not  consistent  to  identify  the  location  of  anything, 
including  a  no-slip  boundary,  to  more  accuracy  than  the  filter  width,  suggesting  that  the  wall  should 
be  filtered  as  well  as  the  turbulence  in  the  interior.  This  approach  provides  a  formal  basis  for  the 
introduction  of  the  wall  stress  in  the  LES  equations,  and  by  including  the  region  exterior  to  the  flow 
domain  in  the  computation,  it  allows  the  wall  stresses  to  be  determined  by  ensuring  that  momen¬ 
tum  and  energy  do  not  “leak”  from  the  flow  domain.  This  simple  and  very  general  modeling  ansatz 
appears  to  be  sufficient  to  provide  LES  boundary  conditions. 

However,  even  with  an  accurate  representation  of  wall  stresses,  a  model  is  needed  for  the  subgrid 
effects  in  the  flow  domain  interior  that  is  valid  near  the  wall  (in  the  log  layer),  where  the  assump¬ 
tions  of  small  scale  isotropy  and  scale  similarity,  on  which  most  subgrid  models  are  based,  are  not 
valid.  Fortunately,  the  optimal  LES  approach  pursued  over  the  years  in  our  group  is  valid  even  in 
the  absence  of  small-scale  isotropy  and  scale  similarity.  Optimal  LES  models  were  used  in  con¬ 
junction  with  the  filtered  boundary /no  leakage  boundary  model,  to  produce  the  encouraging  results 
described  above. 

Application  of  optimal  LES  to  wall  bounded  turbulence  in  conjunction  with  the  flitered  wall/no¬ 
leakage  wall  stress  model  will  be  enabling  for  LES  of  complex  wall-bounded  flows.  The  effective¬ 
ness  of  the  OLES  formalism  in  this  application  was  demonstrated  previously  by  appealing  to  DNS 
statistical  data  for  required  modeling  inputs,  but  to  use  it  in  applications  requires  that  the  need  for 
empirical  data  be  eliminated.  This  can  be  done  using  theory  for  the  multi-point  correlations  of 
turbulence  [26] .  When  small-scale  isotropy  is  valid,  Kolmogorov  inertial  range  scaling  and  mild 
modeling  assumptions  allow  the  required  correlations  to  be  determined,  and  the  resulting  models 
need  to  be  evaluated  and  tested.  For  anisotropic  inhomogeneous  turbulence,  a  generalization  of 
these  correlation  models  is  needed.  An  approach  based  on  the  structure  tensor  representation  of 
anisotropy  was  proposed  earlier.  Here,  we  generalize  and  evaluate  that  approach  for  use  in  wall- 
bounded  turbulence.  Finally,  previous  applications  of  OLES  to  wall  bounded  flows  were  based 
on  a  spectral  representation  of  the  LES  velocity  fields  [39,  5].  To  support  complex  flow  domains, 
the  finite-volume  formulation  of  OLES  needs  to  be  applied  to  wall-bounded  turbulence.  These 
developments  have  been  pursued  in  preparation  for  application  of  OLES  in  complex  wall-bounded 
turbulent  flows. 


3  Supported  Research 


To  pursue  the  objectives  defined  above,  a  number  of  research  activities  were  pursued  under  the 
current  grant.  These  are  described  briefly  below  and  in  more  detail  in  the  following  subsections, 
and  the  referenced  publications. 


1.  Theoretical  OLES  Refinement  &  Testing:  The  most  difficult  statistical  input  to  determine 
theoretically  for  the  OLES  formulation  is  the  three-point  third-order  velocity  correlation 
tensor.  A  recent  development  by  Chang  &  Moser  [9],  provides  such  a  model,  and  so  here, 
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OLES  based  on  this  theoretical  model  has  been  formulated  and  tested.  Further,  an  asymptotic 
analysis  for  small  LES  filter  width  compared  to  the  large  scales  of  turbulence,  produces  a 
particularly  simple  Finite- Volume  OLES  model,  which  can  be  expressed  as  a  finite-volume 
scheme  representing  unresolved  turbulence. 

2.  Refinement  and  Generalization  of  Anisotropy  Models:  An  important  feature  of  near- wall 
turbulence,  is  that  it  is  highly  anisotropic,  and  inhomogeneous.  The  resulting  anisotropy  in 
the  sub-filter  scale  turbulence  needs  to  be  represented,  and  in  the  OLES  framework,  this  is 
done  by  representing  anisotropy  in  the  multi-point  correlations.  We  use  an  approach  for  this 
based  on  the  RANS  structure  tensors  developed  by  Kassinos  &  Reynolds  [14].  The  resulting 
anisotropy  represented  was  tested  for  its  ability  to  reproduce  the  two-point  correlation  in 
channel  flow,  and  it  is  quite  good.  It  is  also  very  complicated,  and  a  potential  alternative, 
more  ad  hoc ,  approach  is  also  discussed. 

3.  Generalization  and  Testing  of  Finite- Volume  Wall-Bounded  OLES:  Finite- volume  for¬ 
mulations  of  OLES  are  much  to  be  preferred  over  spectral  formulations,  that  have  until  now 
been  used  in  wall-bounded  flows.  Wall-bounded  finite-volume  OLES  thus  needs  to  be  for¬ 
mulated  and  tested.  Here  that  is  done  using  statistical  data  from  DNS. 


3.1  Theoretical  OLES  Refinement  &  Testing 

3.1.1  Finite- Volume  Optimal  LES 


In  finite- volume  OLES,  the  LES  state  variables  are  the  velocities  averaged  over  discrete  volumes. 
The  mapping  (filter)  is  equivalent  to  the  application  of  a  top-hat  filter  followed  by  sampling  on  a 
grid.  The  LES  evolution  equations  are  then  determined  from  the  volume  averaged  Navier-Stokes 
equations  given  by  : 

= -Y.T‘ -Y.v‘ +Y.V‘  <d 

s  s  s 

where  u'-  is  the  velocity  averaged  over  the  volume  v,  Vv  is  the  volume  of  v  and  JF®,  V-  and  V®  are 
the  convective  flux,  pressure  force  and  viscous  flux,  respectively,  for  the  surface  s.  The  sums  in  (1) 
are  over  the  faces  of  the  volume  v.  The  quantities  appearing  in  (1)  are  defined  as: 


uv,  = 

v-J„tt,Sx 

(2) 

Ft  = 

J  UiUjU j  fix 

(3) 

vt  = 

J  pn'l  fix 

(4) 

Vt  = 

f  dUi  s 
/  Uf)r  n>SX 

J  S 

(5) 

where  tq  is  the  turbulent  velocity  and  n®  is  the  outward-pointing  unit  normal  to  the  surface  s.  As 
indicated,  the  integrals  are  over  a  volume  v  or  one  of  the  faces  s  bounding  a  volume.  To  distinguish 
the  simulation  quantities  in  an  LES  from  the  filtered  real  turbulence,  the  symbol  will  be  used 
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to  represent  the  LES  variables.  The  goal,  of  course,  is  for  the  dynamics  and  statistics  of  w \  to 
approximate  those  of  ul-  as  closely  as  possible. 

The  evolution  equation  for  will  be  the  same  as  that  for  u"  (1),  with  the  fluxes  replaced  by 
models.  In  the  context  of  OLES,  the  fluxes  Tf,  V,f  and  V,f  are  to  be  modeled  using  stochastic 
estimation.  For  the  current  development,  we  will  consider  the  limit  of  infinite  Reynolds  number  in 
which  the  viscous  flux  V/  is  negligible.  Further,  the  pressure  force  will  be  treated  as  in  Langford 
&  Moser  [20]  (see  section  3.1.3),  so  the  development  in  this  section  will  focus  on  the  convective 
fluxes  These  fluxes  are  to  be  modeled  in  terms  of  linear  and  quadratic  functions  of  the  LES 
state  variables  <•  At  least  quadratic  dependence  is  required  here  because  the  convective  fluxes 
are  themselves  quadratic  in  the  velocity.  The  estimate  for  the  convective  flux  is  thus  of  the  form: 

~  A:(s)  +  Y  vi)wj1  +  Y  vi,v2 )wt-1wl2  (6) 

V\  V\,V2 

where  Zips)  is  a  constant  term,  £tJ(s,  v)  is  the  linear  estimation  kernel,  and  Qijk(s,vi,v2)  is  the 
quadratic  estimation  kernel.  The  range  of  the  sums  in  (6)  can  be  selected  to  be  as  large  or  small  as 
desired,  with  the  expectation  that  a  larger  range  (a  larger  stencil)  will  produce  more  accurate  results, 
though,  as  the  stencil  grows  to  include  more  distant  and  less  well-correlated  data,  a  diminishing 
return  is  expected  [19].  It  was  found  by  ZLM  [42]  that  a  sum  over  4  or  6  volumes  was  sufficient 
to  get  quite  accurate  LES  results.  In  the  LES  performed  here,  a  stencil  is  used  that  sums  over  4 
volumes  (a  1  x  1  x  4  stencil,  see  ZLM  and  section  3.1.3). 

The  minimum  mean  square  error  between  the  ideal  LES  and  the  estimate  (6)  is  attained  when 
[1,19,39]: 

{?!)  =  A(s)  +  Y  vi)  (“?>  +  Ul,  V2)  ( ufulk 2)  (7) 

V\  V\,V2 

(•77«r>  =  A(s)  (wp)  +  Y  Ul)  (ufuj3)  +  Y  Qijk(s,  ui,  u2)  (uful2^3)  (8) 

Vi  Vi,V2 

<xx3<>  =  + 

VI 

+  Y  Ui,  U2)  (upWpMpC)  (9) 

Vl,V2 

These  equations  can  be  solved  for  the  estimation  coefficients  (kernels)  A%,  Cr)  and  Q%3k-  The 
statistical  correlations  appearing  in  the  estimation  equations  (7-9)  are  needed  as  input  for  the  OLES 
procedure. 


3.1.2  Theoretical  Determination  of  Correlations 


The  correlations  required  for  the  OLES  formulation  are  among  the  LES  state  variables  (volume- 
averaged  velocities)  and  the  modeled  quantities  (convective  fluxes).  To  determine  them  theoreti¬ 
cally,  they  must  be  related  to  turbulence  statistical  quantities  for  which  we  have  theories.  This  can 
easily  be  done,  by  writing  the  correlations  appearing  in  (7-9)  as  integrals  of  multi -point  velocity 
correlations.  Three  correlation  tensors  are  needed  to  determine  the  estimation  coefficients.  Models 
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for  these  tensors  must  be  provided  as  input  to  the  OLES  procedure.  To  develop  these  models, 
a  number  of  assumptions  will  be  made.  First  it  is  assumed  that  on  length  scales  required  for 
the  analysis  (i.e.,  several  filter  widths),  the  turbulence  is  both  homogeneous  and  isotropic.  Using 
homogeneity,  the  required  correlations  are  written: 


R.y(rl)  = 

«(x)w'(xl)) 

(10) 

Tijfc(rl,  r2)  = 

(w'(x)w'(xl)i4(x2)) 

(11) 

FjiW(rl,r2,r3)  = 

(u'i  (x)u'-  (xl)u'fc  (x2  )u[  (x3) ) , 

(12) 

with  the  dependence  expressed  in  terms  of  spatial  separations,  ri  =  x  —  x/. 

Second,  it  is  assumed  that  the  spatial  separations  are  small  enough  to  be  within  the  Kolmogorov 
inertial  range  and  that  the  Reynolds  number  based  on  the  filter  width  is  sufficiently  large  that  it  can 
be  considered  to  be  infinite.  With  these  assumptions,  the  Kolmogorov  [16,  15]  expressions  for  the 
second  and  third-order  longitudinal  structure  functions  are  valid: 

S2(r)  =  <(u||(x)  -W||(xl))2)  =Ce2/3r2/3  (13) 

S3(r)  =  ((tt||(x)  -«||(xl))3)  =  -\er  (14) 

where  r  =  |r|  is  the  magnitude  of  the  separation  vector,  u\\  is  the  velocity  component  in  the 
direction  of  the  separation  vector  and  C  is  the  Kolmogorov  constant.  From  these  expressions  and 
isotropy  and  continuity  constraints,  expressions  for  the  two-point  second-  and  third-order  correla¬ 
tions  can  be  derived: 


%(r)  =  u28ij  +  ^e2/3r_4/3(rirj  -  4(r)%)  (15) 

T?:ifc(0,  r)  =  ^  (^ijrk  -  (16) 

The  result  for  the  second-order  correlation  is  well  known.  The  expression  for  the  third-order 
correlation  is  less  common,  but  it  is  a  direct  consequence  of  the  4/5  law  (14)  and  the  general 
isotropic  form  derived  by  von  Karman  &  Howarth  [40]  (see  also  [3,  24]).  Further,  the  analysis 
leading  to  (16)  is  implicit  to  the  derivation  of  the  4/5  law.  The  two-point  third-order  correlation  in 
(16)  is  precisely  the  correlation  needed  to  compute  the  correlation  of  fluxes  with  volume  averaged 
velocities.  However,  the  more  general  third-order  three-point  correlation  is  also  needed  for  the 
correlation  of  three  volume  averaged  velocities  (see  below). 

Finally,  to  determine  an  expression  for  F^,  we  invoke  the  quasi-normal  approximation,  which 
states  that  the  fourth-order  cumulants  are  zero.  This  implies  that  F  can  be  expressed  in  terms  of 
the  two-point  correlation  R: 


Fjjfc/(rl,  r2,  r3)  «  R^(iT)Rw(r3  -  r2)  +  Rifc(r2)R.,7(r3  -  rl)  +  Ra(r3)Rjfc(r2  -  rl)  (17) 

The  quasi-normal  approximation  is  infamous  in  turbulence  because  it  is  well-known  to  result  in 
unrealizable  spectra  when  used  in  two-point  closure  models  [17,  29,  30,  21].  However,  its  appli¬ 
cation  here  in  optimal  FES  modeling  can  result  in  no  such  catastrophe  because  the  FES  equations 
being  solved  are  not  for  statistical  quantities  to  which  realizability  constraints  apply.  Similarly, 
the  modifications  added  to  the  quasi-normal  approximation  to  assure  realizability  in  two-point 
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Figure  1:  Relative  error  <f>ijtki  in  the  quasi-normal  approximation  to  the  non-zero  elements  of  the 
fourth-order  correlation  r,  r),  with  separation  in  the  x\  direction.  To  obtain  0,  the  error  in 

F  is  normalized  by  [F^^F^lO,  r,  r)] 1/2.  Due  to  isotropy  the  2  and  3  indices  can  be  swapped  with 
identical  results. 


closures,  such  as  “eddy  damping”  and  Markovization  in  EDQNM  [21],  are  applied  to  the  evolu¬ 
tion  equation  for  third-order  correlations.  Since  no  such  correlation  evolution  equations  are  being 
solved  here,  these  refinements  are  not  applicable.  Finally,  the  quasi-normal  approximation  is  gen¬ 
erally  quite  accurate  in  isotropic  turbulence  [38].  For  example,  the  error  was  measured  in  a  forced 
isotropic  turbulence  DNS  at  Re A  =  164  [19],  and  is  shown  in  figure  1.  The  small  magnitude 
of  these  errors  does  not  necessarily  imply  that  they  are  dynamically  insignificant.  Ultimately, 
the  performance  of  the  resulting  models  is  the  most  important  measure  of  the  accuracy  of  these 
approximations  (see  section  3.1.3  for  example  LES  results). 

With  the  exception  of  the  third-order  three-point  correlation,  all  the  correlations  required  for  the 
OLES  model  have  now  been  defined.  Unfortunately,  the  simple  theories  employed  here  are  not 
sufficient  to  determine  an  expression  for  the  three-point  third-order  correlation.  Indeed,  just  writing 
down  the  most  general  isotropic  form  satisfying  continuity  constraints  is  difficult  [33].  Based  on 
this  most  general  form,  Chang  &  Moser  [9]  developed  a  model  for  T^r1,  r2)  for  stationary, 
incompressible,  homogeneous,  isotropic  turbulence  for  separations  r1  and  r2  in  the  inertial  range, 
and  under  the  assumptions  of  small  scale  isotropy  and  infinite  Reynolds  number.  This  model  is 
algebraically  very  complex,  with  more  than  758  terms,  so  it  will  not  be  written  here.  Computer 
programs  to  evaluate  the  tensor  numerically  are  available  at  http://turbulence.ices.utexas.edu. 


3.1.3  Testing  Theoretical  Optimal  LES 


To  assess  the  validity  of  the  theoretically  determined  correlations  for  use  in  OLES  models,  several 
of  the  simulations  performed  by  ZLM  [42]  were  repeated  with  the  theoretically  determined  corre¬ 
lations  in  place  of  the  DNS -determined  correlations  used  by  ZLM.  The  details  of  these  simulations 
are  given  briefly  below. 
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The  simulations  were  performed  for  forced  isotropic  turbulence,  using  a  finite-volume  OLES  for¬ 
mulation  like  that  described  here  and,  in  more  detail,  in  ZLM  [42].  As  described  in  section  3.1.1, 
three  fluxes  need  to  be  estimated.  The  convective  flux  T ,  the  pressure  force  V  and  the  viscous  flux 
V.  In  these  simulations,  the  viscous  fluxes  are  zero,  consistent  with  the  infinite-Reynolds  number 
assumption  used  for  the  convective  fluxes.  As  described  by  Langford  &  Moser  [20]  and  as  imple¬ 
mented  by  ZLM  [42],  the  pressure  force  is  determined  by  imposing  an  approximate  divergence- 
free  constraint,  which  is  consistent  with  the  second-order  staggered  grid  finite-volume  divergence 
operator.  This  is  an  optimal  representation  of  divergence  [20],  which  minimizes  the  expected  error 
incurred  by  imposing  an  approximate  divergence-free  constraint. 

The  primary  model  to  be  considered  here  is  that  for  the  convective  fluxes  T .  The  theoretical  models 
described  in  section  3.1.2  for  the  multi-point  correlations  were  integrated  numerically  to  determine 
the  correlations  used  in  the  estimation  equations  (7-9). 

The  integrals  need  to  be  performed  for  each  combination  of  volumes  and  surface  in  the  stencil  and 
for  each  velocity  component  in  the  stencil.  The  stencil  used  here  is  a  1  x  1  x  4  simple  stencil  on 
a  staggered  grid,  which  is  stencil  S4  as  described  by  ZLM  [42].  This  stencil  was  selected  because 
ZLM  found  that  it  is  the  smallest  that  yields  good  a  posteriori  results.  Because  the  stencil  is  defined 
on  a  staggered  grid,  its  definition  is  somewhat  complicated.  See  Appendix  B  in  [42]  and  [26]  for  a 
complete  description. 

Since  the  advent  of  the  dynamic  procedure  in  LES  modeling  [13],  it  has  been  common  practice 
to  seek  ways  in  which  LES  model  parameters  can  be  determined  from  the  LES  in  which  they  are 
being  used.  This  has  the  advantage  of  eliminating  adjustable  model  constants  and  allowing  the 
model  to  respond  to  the  particular  details  of  the  flow.  In  the  OLES  models  developed  here,  only 
those  model  parameters  that  are  clearly  flow-dependent  (i.e.  velocity  variance  u 2  and  dissipation 
rate  e)  are  treated  dynamically,  while  the  other  constants  that  appear  in  the  formulation,  such  as 
the  Kolmogorov  constant,  which  is  given  its  usual  value  C  ~  2,  are  treated  as  universal  constants. 
Setting  these  quantities  directly  reduces  both  the  complexity  of  the  simulations  and  the  chance  that 
the  process  of  determining  the  constants  will  introduce  spurious  dynamics  into  the  simulation. 

Using  the  correlations  derived  from  our  theory,  it  is  straightforward  to  dynamically  determine 
estimates  of  u2  =  2k/ 3  and  e  in  a  running  LES.  That  is  what  is  done  here,  with  the  values  being 
set  at  each  time  step  (see  [26]  for  details)  The  error  in  the  dynamically  determined  value  for  the 
dissipation,  eest ,  computed  using  filtered  data  from  a  DNS  at  Re \  =  164  [19]  is  shown  in  figure  2. 
When  the  separation  r  is  in  the  inertial  range,  eest  is  within  a  percent  of  the  value  determined 
directly  from  the  DNS  (eDNs)- 


3.1.4  Asymptotic  Optimal  LES  Models 


The  estimation  equations  can  be  solved  asymptotically  for  small  7  =  A e/u3.  Use  of  the  low¬ 
est  order  asymptotic  kernels  simplifies  the  LES  model  because  it  removes  the  7  dependence  of 
the  kernels.  Then,  even  when  u2  and  e  are  being  determined  dynamically,  the  scaled  kernels  do 
not  change.  Indeed,  the  lowest  order  quadratic  kernel  Q°  is  scaled  only  with  A,  so,  even  when 
unsealed,  the  kernel  only  depends  on  the  geometry  of  the  cells  being  used  to  estimate  the  fluxes, 
and  not  u2  and  e.  In  this  way,  the  asymptotic  quadratic  kernel  can  be  thought  of  as  a  finite-volume 
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Figure  2:  Relative  error  in  estimating  the  dissipation  from  the  longitudinal  third-order  structure 
function  of  a  finite-volume  filtered  velocity  field.  Filtered  structure  functions  were  computed  by 
filtering  the  DNS  at  Re \  =  164  of  [19].  The  filter  was  defined  on  a  cubical  finite  volume  grid  of 
the  size  noted,  on  a  periodic  domain  of  size  2tt.  Each  finite  volume  is  of  size  A,  and  r  is  the  offset 
between  the  volumes  used  to  compute  the  structure  function. 


Vi 

^nnOA’l)  C^(S,V  l) 

2 

1 

0.05764  0.02623 

-0.39059  -0.12768 

Table  1 :  Values  of  the  elements  of  £°  as  determined  for  the  1x1x4  stencil,  with  volume  labels  as 
defined  in  figure  3.  The  value  of  £°  for  volumes  not  listed  here  are  determined  from  the  symmetry 

scheme  for  the  quadratic  terms  that  is  consistent  with  the  statistics  of  turbulence.  In  addition  to 
exploring  the  performance  of  theoretical  OLES  using  the  full  kernels,  determined  without  asymp¬ 
totic  approximation  (the  linitc-y  kernels),  the  performance  of  the  lowest  order  asymptotic  models 
will  be  evaluated. 

The  asymptotic  solution  for  the  kernels  is  easier  to  understand  when  applied  to  the  “simple”  opti¬ 
mal  models  employed  here.  Because  of  isotropy,  we  need  only  consider  the  velocity  component 
normal  to  the  face  through  which  the  flux  is  being  estimated  and  a  generic  velocity  component 
tangential  to  the  face  being  considered.  In  what  follows,  these  components  will  be  denoted  with 
subscript  n  and  t  respectively  (repeated  n  and  t  do  not  imply  summation).  For  simple  OLES  mod¬ 
els,  flux  of  normal  momentum  is  estimated  in  terms  of  quadratic  products  of  the  form  unun  and 
linear  dependence  on  un  only.  Tangential  fluxes  are  estimated  with  quadratic  and  linear  terms  of 
the  form  unut  and  ut,  respectively.  As  a  consequence,  the  only  quadratic  kernel  elements  are  Qnnn 
and  Qtnt,  and  the  linear  kernel  elements  are  Caa,  where  the  subscript  a  can  be  either  n  or  t,  and 
no  summation  is  implied. 

Values  of  the  scaled  asymptotic  kernel  elements  are  given  in  table  1  and  2  for  the  1x1x4  stencil. 
In  these  tables,  the  volume  labels  are  as  defined  in  figure  3,  for  estimates  of  the  flux  through  the 
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Figure  3:  Volumes  (and  their  numbering)  in  a  1  x  1  x  4  staggered  grid  stencil  to  estimate  the  flux 
through  the  surface  between  volumes  -1  and  1.  For  the  normal-component  flux  (n),  the  volumes 
with  solid  outlines  are  on  the  n-componcnt  mesh,  and  the  dashed-outline  volumes  are  not  used. 
For  the  tangential-component  flux  (£),  the  solid  volumes  are  on  the  /-component  mesh,  and  the 
dashed  volumes  are  on  the  n-component  mesh.  In  tables  1,  2,  the  volumes  are  referred  to  by  the 
numbers  shown  here,  but  no  distinction  is  made  between  0+  and  CT  because  the  values  associated 
with  these  volumes  are  the  same. 


V\ 

V2 

Q° 

{s,Vi,V  2) 

2 

2 

0.00913 

2 

1 

-0.11485 

2 

-1 

-0.21432 

2 

-2 

0.04386 

1 

1 

0.21467 

1 

-1 

1.16690 

Vl 

V2 

Qtnt 

s,v  1,V2) 

0 

2 

-0.06617 

0 

1 

0.31617 

Table  2:  Values  of  the  elements  of  Q°  as  determined  for  the  1x1x4  stencil,  with  volume  labels  as 
defined  in  figure  3.  The  value  of  £°  for  volumes  not  listed  here  are  determined  from  the  symmetry 

Q°ana(S,~V  1,  =  Q°a  na{s,V1,V2). 
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Figure  4:  Three-dimensional  energy  spectra  (a)  and  third-order  structure  functions  (b)  from  OLES 
of  isotropic  turbulence  at  infinite  Reynolds  number  using  the  finite-7  kernels,  with  resolutions 
ranging  from  163  to  1283  (7  «  0.17  to  7  ps  0.02,  respectively).  The  solid  lines  in  both  plots  are 
determined  from  Kolmogorov  theory.  In  (a),  the  two  solid  lines  are  a  k~ 5//3  slope  (shallow),  and 
the  result  of  filtering  a  k~ 5/3  spectrum.  In  (b)  the  straight  line  is  S3  =  —  |er,  and  the  other  solid 
line  is  the  structure  function  of  the  filtered  velocity  determined  from  I3. 


cell  face  between  volumes  -1  and  1.  The  results  of  performing  an  LES  with  the  asymptotic  kernels 
are  compared  to  that  for  the  finite-7  kernels  in  section  3.1.5. 


3.1.5  LES  Results 


The  performance  of  the  theory-based  OLES  models  developed  here  is  evaluated  in  isotropic  tur¬ 
bulence  at  infinite  Reynolds  number,  simulated  in  a  periodic  cube  with  sides  of  length  L  =  2tt. 
The  three  dimensional  energy  spectra  and  third-order  structure  function  obtained  using  a  1  x  1  x  4 
stencil  [42]  on  a  staggered  grid  are  shown  in  figure  4  for  different  grid  resolutions  (163  to  1283)  cor¬ 
responding  to  finite-volume  filter  widths  varying  from  A  p=i  0.39  to  0.05.  In  these  forced  isotropic 
turbulence  simulations,  the  energy  is  injected  at  the  rate  e  ~  62.3468,  and  the  simulations  result  in 
a  vS  ~  28,  yielding  7  ss  0.17  to  0.02  for  the  163  to  1283  grid  sizes,  respectively.  In  all  the  spec¬ 
tra  presented  here,  the  wavenumber  k  is  normalized  by  kmin  =  2n/L,  and  the  three  dimensional 
energy  spectrum  E(k)  is  normalized  by  e2/3/c~ f/3. 


For  all  filter  widths,  the  high  wave-number  portion  of  the  LES  energy  spectrum  exhibits  a  slope 
consistent  with  the  filtering  of  a  k~ 5/3  inertial  range.  For  larger  grid  sizes  and  at  lower  wavenum¬ 
bers,  the  shallower  k~ 5//3  slope  is  evident,  and  this  range  becomes  longer  with  increased  resolution. 
This  is  the  converged,  resolution-independent  part  of  the  solution.  The  shifting  of  the  filtered  part 
of  the  spectrum  to  higher  wavenumbers  with  increasing  resolution  is  expected,  since  the  filter  width 
is  being  decreased.  Similarly,  the  third-order  longitudinal  structure  function  agrees  very  closely 
with  the  theoretical  filtered  structure  function  obtained  by  integrating  the  third-order  three-point 
correlation  model  [9]  over  the  finite  volumes.  As  expected,  as  the  resolution  increases,  the  sim¬ 
ulation  structure  functions  approach  the  —  |er  dependence  expected  in  the  inertial  range.  This 
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Figure  5:  Three-dimensional  energy  spectra  (a)  and  third-order  structure  functions  (b)  from  OLES 
of  isotropic  turbulence  at  infinite  Reynolds  number  using  both  finite  7  and  asymptotic  kernels 
(signified  with  “A”),  with  resolution  32s  and  1283  (7  =  0.08  and  7  =  0.02  respectively).  Solid 
lines  in  (a)  and  (b)  are  as  in  figure  4. 


suggests  that  LES  is  correctly  representing  the  inertial-range  energy  transfer.  In  short,  the  behavior 
exhibited  in  figure  4  is  precisely  what  is  expected,  indicating  that  the  theory-based  OLES  model 
is  representing  the  effects  of  the  unresolved  scales  at  infinite  Reynolds  number  in  a  consistent, 
resolution  independent  way. 

The  LES  results  presented  in  figure  4  were  obtained  using  the  finite-7  OLES  kernels,  but  as  indi¬ 
cated  in  section  3.1.4,  it  will  generally  be  more  convenient  to  use  the  small-7  asymptotic  kernels. 
A  comparison  of  LES  results  from  the  finite-7  and  asymptotic  kernels  is  shown  in  figure  5.  At 
both  1283  and  323  resolutions,  corresponding  to  7  «  0.02  and  0.08,  respectively,  the  spectra  and 
structure  functions  from  the  two  cases  are  indistinguishable. 

A  finite  Reynolds  number  case  was  simulated  using  OLES  to  gauge  the  performance  of  the  optimal 
models  compared  to  a  standard  model  (dynamic  Smagorinsky).  Spectra  and  structure  functions 
from  these  simulations  are  shown  in  figure  6  along  with  that  from  filtered  DNS.  Results  for  both 
1x1x2  and  1x1x4  stencils  on  a  staggered  32s  grid  are  shown.  With  the  dynamic  Smagorinsky 
model,  the  two  stencils  correspond  to  second-  and  fourth-order  schemes.  Comparing  the  1x1x4 
stencil  optimal  model  to  the  corresponding  dynamic  model,  a  sharp  roll-off  of  the  spectrum  at  high 
wavenumbers  is  evident  with  the  dynamic  model,  but  not  present  in  the  optimal  model.  The  optimal 
model  apparently  yields  a  better  treatment  of  the  model  dissipation  at  high  wavenumbers.  However, 
the  1  x  1  x  2  stencil  optimal  model  does  have  the  sharp  roll-off.  The  optimal  model  spectra  also 
have  a  somewhat  steeper  slope  than  either  the  filtered  DNS  or  the  dynamic  model,  but  this  slope 
is  consistent  with  the  filtered  k~5^3  spectrum.  Optimal  models  for  these  stencils  that  are  based  on 
the  DNS  correlations  do  not  exhibit  this  difference  in  slope  of  the  spectrum  between  the  model  and 
filtered  DNS  (see  ZLM  [42]).  It  appears  that  the  optimal  model  results  reflect  the  infinite-Reynolds 
number  theory  on  which  the  correlations  are  based,  which  is  not  strictly  applicable  to  this  moderate 
Reynolds  number  case.  For  the  third-order  structure  function,  the  optimal  models  are  somewhat 
more  accurate  in  the  mid-range  of  separations  than  the  dynamic  Smagorinsky.  Otherwise,  the 
model  results  are  comparable. 
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Figure  6:  Three-dimensional  energy  spectra  (a)  and  third-order  structure  function  (b)  from  LES 
of  isotropic  turbulence  on  a  323  finite-volume  grid  at  Re A  =  164  using  the  dynamic  Smagorinsky 
model  and  finite-7  optimal  models,  with  1x1x2  and  1x1x4  stencils.  Also  shown  is  a  filtered 
DNS.  Numbers  on  the  curve  labels  indicate  the  stencil  size,  2  signifies  the  1x1x2  stencil,  4  is 
the  1x1x4  stencil. 


At  high  Reynolds  numbers  these  optimal  models  should  yield  accurate  results  (as  seen  in  figure  4), 
and  will  generally  be  superior  to  the  dynamic  model  due  to  the  elimination  of  the  sharp  spectral 
roll-off  near  the  cutoff. 


3.1.6  Role  of  the  Three-Point  Third-Order  Correlation 


The  most  difficult  modeling  task  involved  in  formulating  the  models  presented  here  is  the  represen¬ 
tation  of  the  three-point  third-order  correlation.  The  representation  for  this  quantity  was  feasible 
[9]  only  because  of  the  assumption  of  small-scale  isotropy.  In  more  complex  situations  in  which 
the  small-scale  isotropy  assumption  would  not  be  valid,  detailed  modeling  of  the  three-point  third- 
order  correlation  will  be  much  more  difficult,  if  not  infeasible.  It  is  interesting,  then,  to  consider 
how  important  this  quantity  is  to  the  modeling  and  thus  how  well  it  needs  to  be  known. 

The  primary  use  of  the  three-point  correlations  is  in  the  terms  that  couple  the  equations  for  the 
quadratic  and  linear  kernels.  The  coupling  is  most  important  in  the  equation  for  the  linear  kernels, 
as  the  coupling  term  in  the  quadratic  equation  is  second  order  in  y2/3.  The  equation  for  the  linear 
kernels  can  be  interpreted  as  a  condition  requiring  that  the  contributions  of  the  model  to  the  filtered 
two-point  correlation  be  the  same  as  that  for  true  turbulence,  at  least  for  separations  included  in 
the  kernel  stencil  (see  [19]).  The  critical  role  of  J3,  the  integral  of  the  three-point  correlation,  then, 
is  to  measure  this  contribution  of  the  quadratic  part  of  the  model.  In  particular,  it  measures  the 
contribution  of  the  quadratic  model  term  to  the  transfer  of  energy  to  the  small  scales.  However, 
the  contribution  to  the  dissipation  (transfer  to  small  scales)  of  the  asymptotic  quadratic  model  term 
is  actually  quite  small  (about  2%,  see  [26]).  If  this  were  generally  indicative  of  the  contribution 
of  the  I3  terms  to  the  estimation  equations,  then  it  would  be  reasonable  to  neglect  them,  greatly 
simplifying  the  modeling  problem  by  making  it  largely  unnecessary  to  develop  a  model  for  the 
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Figure  7:  Three-dimensional  energy  spectra  (a)  and  third-order  structure  functions  (b)  from  OLES 
of  isotropic  turbulence  at  infinite  Reynolds  with  resolution  of  323  (7  =  0.08).  Compared  are  results 
using  the  finite  7  kernels,  asymptotic  kernels  (signified  with  “A”)  and  asymptotic  kernels  generated 
by  neglecting  the  / 3  terms  in  the  estimation  equations  (signified  with  ”N”).  Solid  lines  in  (a)  and 
(b)  are  as  in  figure  4. 


three-point  third-order  correlation. 

To  explore  this  possibility,  a  set  of  asymptotic  kernels  were  determined  by  neglecting  the  I3  terms 
in  the  asymptotic  estimation  equations.  Large  eddy  simulations  with  these  modified  kernels  were 
performed,  and  the  results  for  the  three-dimensional  spectrum  and  third-order  structure  function 
are  shown  in  figure  7.  The  spectrum  in  the  simulation  that  neglected  the  I3  terms  differs  from  that 
obtained  with  the  usual  asymptotic  kernels  in  the  highest  octave  of  wavenumbers,  with  an  up-turn 
followed  by  a  sharper  roll-off.  This  demonstrates  the  role  of  the  I3  terms  in  shaping  the  scale- 
distribution  of  the  dissipation  provided  by  the  linear  part  of  the  model,  in  addition  to  the  overall 
dissipation  rate  described  above.  The  effect  of  neglecting  the  J3  terms  on  the  simulated  third-order 
structure  function  is  pronounced  only  at  separations  greater  than  4 A,  for  reasons  that  are  not  clear. 

The  importance  of  these  differences  to  a  practical  LES  can  be  debated.  In  general,  whether  or  not 
the  additional  modeling  error  introduced  by  neglecting  the  I3  terms  is  acceptable  will  depend  on 
the  goals  of  a  particular  simulation.  The  modest  impact  on  the  LES  results  in  this  case  suggests  that 
neglecting  the  I3  terms  or  crudely  modeling  the  underlying  three-point  third-order  correlation  (e.g. 
to  account  for  anisotropy)  may  be  a  viable  strategy  in  more  complex  flow  situations.  However, 
one  must  be  cautious  in  this  conclusion  since  in  more  complex  situations,  the  quadratic  operator 
may  not  behave  as  well,  requiring  a  more  tailored  correction  with  the  linear  term,  which  would  be 
determined  through  the  I3  terms. 


14 


3.2  Refinement  and  Generalization  of  Anisotropy  Models 

For  a  given  velocity  field  u(x),  the  two-point  correlation  tensor  F!tJ  is  given  by 

#ij(x,  r)  =  (m'(x)m'-(x  +  r))  ,  (18) 

where  u'  =  u  —  (u)  is  the  fluctuating  velocity  field.  Two-point  correlations  of  velocity  fluctuations 
are  an  important  ingredient  of  the  optimal  LES  formulations  pursued  here,  as  discussed  in  the 
previous  section.  Indeed  a  correlation  representation  with  a  finite  number  of  parameters  can  be 
filtered  and  then  fit  to  the  correlations  of  an  LES.  The  resulting  parameters  then  yield  a  model  for 
the  underlying  unfiltered  correlations.  This  ability  to  infer  the  correlations  of  u'  from  the  statistics 
of  the  filtered  fields  is  fundamental  to  LES,  it  forms  the  basis  of  subgrid  modeling,  and  is  necessary 
if  one  is  to  use  the  results  of  LES  for  analysis  of  turbulent  flows.  It  is  this  need  to  represent  velocity 
correlations  in  turbulent  flows  with  complex  geometries  that  motivates  the  current  work. 

At  any  given  x  location,  one  of  the  key  properties  of  Riq  that  needs  to  be  represented  is  its 
anisotropy,  both  with  regard  to  how  different  H,  j  (x.  0)  is  from  i?^(x,  0)^/3,  (componental 
anisotropy)  and  with  regard  to  the  elongation  of  the  isocontours  of  Raa(x,  r)  for  different  direc¬ 
tions  in  separation  (directional  anisotropy).  Perhaps  the  most  notable  effort  to  represent  this 
anisotropy  was  made  in  the  work  by  Arad  et  al[ 2],  which  we  will  briefly  describe,  before  dis¬ 
cussing  our  motivation  to  pursue  a  different  approach.  In  the  approach  taken  by  Arad  et  al[ 2],  a 
complete  basis  for  Rij  was  formulated  in  terms  of  subspaces  that  are  invariant  to  rigid  rotation  of 
the  frame  of  reference,  i.e.  an  SO(3)  decomposition  was  proposed  for  /?,  ;(r  ).  Thus,  a  given  second 
rank  tensor  function  of  r,  say  Tt](r),  can  be  expanded  as  T^(r)  =  Tt](0)  +  J2qirn  aqim(r)Bfjm(r /r) 
(no  summation  implied  on  repeated  indices).  Here,  B^m(r/r)  are  basis  tensors  depending  only 
on  the  angular  variation  in  r,  while  aqim(r )  are  scalar  functions.  For  a  given  (7,m),  the  set 
{ B'fi'"(r/r):  q  =  1, ...  Nq(l)}  is  generated  from  Y)m(r/r),  the  spherical  harmonics,  and  this  is  a 
finite-dimensional  subspace  (i.e.  Nq(l)  is  finite  for  all  /)  that  is  invariant  to  rotational  transfor¬ 
mations.  In  the  context  of  turbulent  flows,  for  expansions  of  Rij{ r),  a  power-law  form  was  then 
proposed  for  the  scalar  functions,  i.e.  aqim(r )  =  cqimr^l\  and  it  was  speculated  that  if  /i  >  l2  then 
('(/,  )  >  ('(/a).  The  existence  of  a  hierarchy  for  the  power-law  exponents  would  of  course  be  very 
valuable,  because  at  small  r,  the  modes  with  lower  l  would  dominate,  and  therefore  R,:j(i')  could 
be  well  approximated  with  a  truncated  series,  yielding  a  finite-dimensional  representation. 

A  power  law  hierarchy  has  indeed  been  found  in  correlations  constructed  from  DNS  of  homoge¬ 
nous  turbulent  flows  [7].  However,  DNS  data  of  wall-bounded  flows  at  high  Reynolds  number  do 
not  yield  clear  evidence  of  this  hierarchy  amongst  the  anisotropic  (/  >  0)  modes[8].  A  similar 
result  was  obtained  when  we  performed  SO(3)  decompositions  of  R,j  (r)  computed  from  DNS  at 
Rer  =  940  (not  shown  here).  Inhomogenie ty  in  the  wall-normal  (y)  direction  is  a  possible  reason 
for  the  lack  of  hierarchy  of  power-law  exponents  in  wall-bounded  flows.  Kurien  and  Sreenivasan 
[18]  accounted  for  the  r  ~  y  range  of  Rij(r)  while  measuring  exponents  in  a  turbulent  bound¬ 
ary  layer,  and  found  that  £(2)  was  much  larger  than  ('(()),  consistent  with  an  approach  to  isotropy 
at  small  scales.  However,  no  comparison  was  made  amongst  the  l  0  exponents.  Without  a 
hierarchy  in  the  power-law  exponents,  it  is  difficult  to  choose  one  set  of  SO(3)  modes  over  others. 

The  abovementioned  issues  with  the  SO(3)  decomposition  in  the  context  of  wall-bounded  flows 
motivated  us  to  investigate  a  different  approach.  Here,  Rij(r )  is  represented  approximately  in 
terms  of  single-point  tensors,  the  structure  tensors  proposed  by  Kassinos  et  al  [14]  (referred  to 
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below  as  KRR).  These  second-order  tensors  are  single-point  moments  of  derivatives  of  fluctuat¬ 
ing  stream  functions,  which  can  therefore  be  related  to  integrals  of  two-point  correlations  over  r. 
Though  these  tensors  are  single-point  moments,  they  contain  information  about  the  “structure”  of 
the  anisotropy,  e.g.  the  distribution  of  energy  in  different  components  of  velocity,  the  dimension  of 
turbulence,  etc  (discussed  in  [6]).  Also,  in  the  context  of  constructing  a  closure  model  for  Reynolds 
stresses,  it  was  shown  by  KRR  that  it  is  vital  to  represent  the  pressure-strain  correlation  in  terms 
of  structure  tensors  in  order  to  obtain  the  correct  evolution  of  the  Reynolds  stress  components  in 
the  rapid  distortion  limit.  The  experience  of  KRR  with  the  structure  tensors  suggests  that  they 
encode  anisotropy  information  that  is  important  to  the  evolution  of  the  turbulence.  This  motivates 
the  proposed  ansatz  that  the  anisotropy  in  Rij  be  expressed  in  terms  of  structure  tensors  through 
the  use  of  the  theory  of  invariants  [35]  to  construct  the  most  general  linear  form. 

The  obvious  advantage  of  this  approach  over  the  one  by  Arad  et  al[ 2]  is  that  even  the  most  general 
basis  of  tensor  functions  obtained  here  is  finite-dimensional,  so  no  power  law  hierarchy  is  needed 
to  truncate  the  representation.  However,  unlike  the  SO(3)  representation,  the  linear  representaion 
in  terms  of  structure  tensors  is  not  a  complete  basis  for  R,tJ(r).  The  quality  of  the  representation 
depends  on  the  validity  of  the  modeling  ansatz  on  which  it  is  based.  To  evaluate  the  capabilities  and 
shortcomings  of  our  representation,  we  obtain  a  model  for  Rt:  by  fitting  it  to  correlations  obtained 
from  DNS,  and  compare  the  two  correlations  (Sec.  3.2.3). 

The  results  of  this  anisotropy  development  are  described  briefly  below,  and  in  more  detail  in  [6]. 


3.2.1  Background 


Structure  tensors  are  single-point  moments  of  derivatives  of  stream  functions.  The  set  of  struc¬ 
ture  tensors  that  are  nonzero  for  homogeneous  turbulence  are  given  by  the  Reynolds  stress  (or 
componentality)  BtJ,  dimensionality  1),,,  circulicity  Ft]  and  stropholysis  Q*jk.  In  addition  for 
inhomogeneous  turbulence  the  inhomogeneity  tensor  is  non-zero.  See  [14]  for  definitions. 
These  tensors  measure  different  characteristics  of  the  anisotropy,  particularly  the  anisotropy  of  the 
velocity  components  ( B ),  the  anisotropy  of  correlation  lengths  (I)),  the  anisotropy  of  the  vorticity 
components  (F),  the  breaking  of  planar  reflection  symetries  (0)  and  the  anisotropy  arising  due 
to  inhomogeneity  (C).  While  these  structure  tensors  are  single -point  statistical  quantities,  they 
can  be  expressed  in  terms  of  the  two-point  second-order  velocity  correlation,  and  can  therefore  be 
interpreted  as  measures  of  the  anisotropy  of  this  correlation  [6]. 


3.2.2  A  Model  Form  for  the  Homogeneous  Anisotropic  Two-Point  Correlation 


In  formulating  a  structure-tensor  based  anisotropic  model  of  the  two-point  correlation  tensor,  it 
will  be  convenient  to  consider  an  infinite  Reynolds  number  representation  for  R^,  with  a  finite 
Reynolds  number  (viscous)  correction.  In  the  limit  of  infinite  Reynolds  number,  the  inertial  range 
variation  of  the  correlations  (assumed  to  follow  a  power-law)  extends  all  the  way  to  zero  separation, 
where  the  derivatives  of  the  correlation  are  then  discontinuous.  At  finite  Reynolds  number,  the 
discontinuities  at  r  =  0  will  be  “healed”  in  a  region  of  r  with  size  of  order  the  Kolmogorov  scale 
rj,  so  the  finite  Re  correlation  must  diverge  from  the  infinite  Re  correlation  for  small  r  (see  Fig.  8). 
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Figure  8:  Schematic  showing  the  infinite  and  finite-Re  stream- wise  correlations  superimposed  on 
each  other.  r)  coincides  with  Ru( r)  between  re  <  r  <  tl,  hut  due  to  viscous  healing  peels 
off  below  r  =  re.  The  relative  size  of  re  has  been  exaggerated  for  the  sake  of  clarity. 


To  distinguish  the  finite  and  infinite  Reynolds  number  quantities,  a  superscript  v  will  indicate  the 
finite  Reynolds  version,  so  that: 

ity(r)  =  lim  fiTj(r),  (19) 

tie — >oo 

Bn  =  lim  B”  (20) 

Re — >oo  J 

The  inertial  range  over  which  the  infinite  Reynolds  number  model  for  i?^(r)  will  be  valid  extends 
over  re  <  r  <  rL,  where  re  ~  rj  <C  rL  and  rL  <C  £  (£  being  the  length  scale  of  the  flow  geometry), 
and  the  viscous  healing  occurs  in  R^( r)  for  r  <  re,  as  shown  in  Fig.  8.  The  representation  of  R,j 
in  terms  of  structure  tensors  models  the  infinite-Re  Rrj  (r) .  A  weakly  anisotropic  finite  Reynolds 
number  correction  is  then  introduced  for  r  ~  r/.  The  correction  will  be  important  when  we  fit  our 
representation  to  finite-Reynolds  number  correlations  computed  from  DNS. 


In  [6],  it  was  shown  that  the  structure  tensors  encode  important  information  about  the  anisotropy 
of  Rij.  This  information  is  now  assumed  to  be  sufficient  to  reproduce  the  anisotropic  features  of 
the  correlation  in  anisotropic  turbulent  flows.  The  correlation  is  thus  assumed  to  depend  on  the 
structure  tensors  and  the  separation  vector  r.  Since  we  are  striving  to  represent  anisotropy,  it  will 
be  convenient  to  re-express  the  dependencies  in  terms  of  the  anisotropy  tensors  associated  with  B, 
D  and  Q*,  which  are  given  by 


bij 


dij 


Rij  dij 

~(f  ~  T 

Rij  dij 
~q^~  ~3 

Qijk 


(21) 

(22) 

(23) 
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where  q 2  =  Bn.  A  model  for  Ri:/(r  )  of  the  form 

Rij{ r)  =  Bi:j  +  Fij( (r,  b,  d,  Q*)  (24) 

is  thus  sought,  which  is  formulated  so  that  FtJ  goes  to  zero  at  r  =  0.  Furthermore,  the  function  Ft/ 
will  be  required  to  satisfy  the  following  restrictions: 


•  The  model  for  Ri3  is  meant  for  inertial  range  separations  re  <  r  <  rL,  and  our  modeling 
assumption  is  that  the  effects  of  the  large  scale  will  be  negligible  within  this  range.  Any  fit¬ 
ting  to  DNS  data  will  be  done  over  only  this  range.  Therefore,  Fl}  does  not  depend  explicitly 
on  rL,  and  we  do  not  explicitly  model  the  correlation  for  r  >  rL. 

•  The  dependence  on  bij,dij  and  Q*-k  is  linear.  This  assumption  is  consistent  with  the  exact 
linear  relationship  between  the  structure  tensors  and  Rij. 

•  Fij( r,  b.  d.  Q*)  is  invariant  to  proper  rotation  of  the  reference  frame  as  well  as  to  changes  in 
the  handedness  of  the  axes  of  the  reference  frame. 

•  The  symmetry  Rij(r)  =  Rji(— r)  is  satisfied,  which  is  exactly  true  for  homogeneous  turbu¬ 
lence. 


The  most  general  linear  representation  Fi3( r,  b.  d.  Q*)  is  constructed  using  the  invariant  theory  of 


tensors  [35]  to  obtain  the  following  form  for  R,t] : 

Rij( r)  =  Btj  +  R*j{ r)  +  RbtJ{ r,  b)  +  R^( r,  d)  +  i?g(r,  Q*),  (25) 

where, 

R{j(?)  =  h{r)8ij  +  f2{r)rirj>  (26) 

Rbij(r ,  b)  =  f:>,(r)b;j  +  [/4(r)%  +  f:,(r) v  ■  b  ■  r  (27) 

+/e(^)[A( r  •  h)j  +  r^v  •  b);], 

Rij(r,  d)  =  f7(r )dij  +  [f8(r)Si:j  +  /9(r)ryr.,]r  •  d  •  r  (28) 

+/io('r)['ri(r  •  d)j  +  r j{r  ■  d)j], 

R?j(r,  Ql  =  fn(r)hmkQ*klj  +  ejmkQluhrm  (29) 

+fi2(r)[rjeink  +  riejnk\Q*klrnrirmrn. 


Here  j\  (r) — f  12(f)  are  scalar  functions  of  r.  The  number  of  free  scalar  functions  f,(r)  is  reduced 
by  invoking  constraints  imposed  by  continuity  and  by  self-consistency  of  the  representation.  In 
particular,  a  constraint  is  imposed  that  terms  involving  d,  for  example,  not  contribute  to  the  com- 
ponentality  B  or  stropholysis  Q  of  the  model  correlation,  and  similarly  for  the  other  terms. 


The  resulting  constraints  are: 


fi  +  r2f2  +  4  rf2 

=  0, 

(30) 

fi  +  6  rh  +  r2/s  +  /e 

=  0, 

(31) 

2  r/4  +  /s  +  5  rfe  +  r2/e 

=  0, 

(32) 

fs  +  6  r/9  +  r2/'  +  f[0 

=  0, 

(33) 

2  r/8  +  /y  +  5r/i0  +  r2f[0 

=  0, 

(34) 

f2  f'\2  +  2'r  fl2  +  fu 

=  0, 

(35) 
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for  continuity,  and: 


2b- ■ 

-^(15g  +  9rg'  +  r2g"  +  35r2f5  +  13r3f'  +  r4f£)  =  0,  (36) 

(30  f'/r  +  15/y  +  40r2/g  +  20r3/'  +  2  r4/'' 

+60/10  +  60r/(0  +  10r2/")  =  0,  (37) 

for  self-consistency,  where  g  =  3/4  +  2 /6  (see  [6]  for  details). 

Because  of  the  linearity  of  the  constraints  and  the  assumed  linear  dependence  of  Rt-j  on  the  struc¬ 
ture  tensors,  the  continuity  and  self-consistency  constraints  described  above  act  individually  on 
R1 ,  Rb,  Rd  and  HP,  without  coupling  between  them.  For  each  of  these  terms  in  the  anisotropy 
model,  the  constraints  leave  one  scalar  function  undetermined.  We  assume  here  that  these  func¬ 
tions  are  power-laws  in  r,  consistent  with  the  expected  functional  form  of  the  correlation  in  the 
inertial  range.  The  constraints  then  require  that  each  of  the  four  terms  in  the  model  for  Rij  has 
an  overall  power-law  dependence  on  the  separation  magnitude  r,  with  a  single  exponent  for  each 
term.  Call  these  four  power-law  exponents  pr,  pb,  Pd  and  pQ.  The  scalar  functions  fjr)  then  have 
the  form  /*(r)  =  alPp,'^z',\  where  a  is  one  of  I,  b,  d,  Q,  depending  on  the  term  in  which  f  ,  appears 
and  Zi  is  a  positive  integer,  which  is  the  net  power  of  r  multiplying  f,  in  (26)-(29).  For  example, 
£3  =  0,  Z4  =  2  and  z5  =  4.  The  coefficients  a,:  are  also  constrained,  so  that  only  four  of  them  are 
independent.  We  choose  to  specify  a  1,  a3,  a7  and  an ■  Given  a  set  of  power-law  exponents  pa,  the 
ratios  of  the  rest  of  the  coefficients  to  these  four  coefficients  are  fixed  by  the  constraints.  Thus  we 
can  rewrite  Eqns.  (26)-(29)  as: 

Rij(r)  =  Bi:j  +  ARi:j(r),  (38) 

where 

Ai?y(r)  =  apr^-2  r2c%  +  —  ryry 

CL\ 

+a3rPb~4\r4bij  +  [— 4,-r2  +  — ryryjr  •  b  •  r  +  — r2[r*(r  •  b)j  +  ry(r  •  b)j]| 
l  a3  a3  a3  J 

+a7rPd~4(r4:dij  +  [— bif2  +  —  rpr^r  ■  d  ■  r  +  —  r2[ri(r  ■  d  )j  +  rj(  r  ■  d)j]| 

l  Oj7  CL  7  fl7  J 

+anrPQ~4{\eirnkQ*klj  +  ejmkQ*kH]r2rirm 

H  —  Ij" j^-ink  T  ’C’iCjnk]  Q klmP l  P |  •  (39) 

an  J 

The  above  relation  implies  that  if  we  fix  (B,  b,  d,  Q*}  then  there  are  8  degrees  of  freedom  (dof) 
in  the  representation,  given  by  the  four  pa  and  ai,  a3,  a7,  au. 

At  zero  separation,  the  infinite  Reynolds  number  inertial-range  model  in  Eqn.  (38)  is  equal  to 
the  infinite  Reynolds  number  Reynolds  stress  Bij.  While  the  finite  and  infinite  Reynolds  number 
correlations  will  match  in  the  inertial  range,  the  Reynolds  stress  TC  will  not,  as  indicated  in  Fig.  8. 
It  is  important  to  account  for  this  difference  A B^  =  B,,  —  /j'(  when  fitting  the  inertial  range 
correlation  model  to  data  from  finite  Reynolds  number  turbulence  as  is  done  in  Sec.  3.2.3,  or  when 
using  the  model  to  evaluate  the  Reynolds  stress  in  a  turbulent  flow  (after  fitting  to  an  LES  for 
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example).  While  ||A5y||  is  usually  small  compared  to  Bkk  at  large  Re,  an  incorrect  estimate  of 
Bij  (for  instance,  using  13,, j  =  B0  shifts  the  correlation  by  A B,tj  over  all  separations,  including 
the  inertial  range.  Thus,  a  good  estimate  of  A B^  is  needed  even  to  match  the  model  to  finite 
Reynolds  number  correlations  in  the  inertial  range. 

As  discussed  above,  Rij{ r)  and  ///(T)  will  correspond  in  the  inertial  range,  but  for  sufficiently 
small  separations  (r  <  re),  viscosity  “heals”  the  derivative  discontinuity  that  would  be  present  if 
the  inertial  range  behavior  extended  to  r  =  0.  In  this  viscous  region,  the  finite  Re  correlation  can 
be  represented  as  a  second  order  Taylor  series  around  zero,  resulting  in  a  composite  representation: 

™  +  MijMrkri  for  r  <  re 

'  [  Rij(r)  for  r  >  re  ’ 

where  is  the  tensor  Taylor  series  coefficient.  For  simplicity,  the  matching  will  be  done  in 
two  steps.  First,  re  will  be  determined  by  matching  the  isotropic  parts  of  the  inner  and  outer 
approximations,  assuming  that  the  isotropic  parts  of  both  approximations  dominate  at  these  small 
separations.  For  this  to  be  a  self-consistent  assumption,  any  anisotropy  in  the  dissipation  rate 
tensor  el3  must  be  weak,  and  the  anisotropic  power  law  exponents  {pb,Pd,PQ}  must  be  greater 
than  pr.  Furthermore,  it  is  assumed  that  the  isotropic  part  of  the  inertial  range  correlation  satisfies 
Kolmogorov’s  2/3  law,  that  is  pr  =  2/3  and  a\  =  — 2Cfce2//3/3,  where  e  =  e^/2  is  the  rate  of 
dissipation  of  kinetic  energy.  Equating  the  isotropic  parts  of  the  viscous  and  inertial-range  models 
at  r  =  re,  one  obtains: 
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The  value  of  re  is  determined  by  matching  the  coefficients  of  (0  obtain: 


=  (5Ck)~f'ri. 


(41) 


Where  q  =  (u3 /e)1/4  is  the  Kolmogorov  length  scale  and  Ck  is  the  Kolmogorov  constant.  Consis¬ 
tent  with  experimental  measurements  [36]  a  constant  value  of  Ck  =  2  is  assumed  here. 


In  the  second  step  in  the  matching  process,  the  correction  to  BtJ  is  determined  by  matching  the 
inner  and  outer  models  at  re  in  an  integral  sense: 

J  [B^  +  ARij(r)}  | r=rdZ  =  j  [B^  +  Mijklrkri\  \ r=rdZ.  (42) 

Making  use  of  the  fact  that 


d^R^v) 
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and  substituting  Eqn.  (43)  and  Eqn.  (39)  into  Eqn.  (42)  yields: 


B{j  —  Bi3  +  A  B^, 


(43) 


(44) 
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Table  3:  Specifics  of  the  turbulent  channel  flow  field  (dimensionless) 


ReT  940 

Resolution  3072  x  385  x  2304 

Box  Size  (lx/h  x  lz/h)  8n  x  3n 


where 
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Here  pr  =  2/3  is  now  a  fixed  parameter  in  accordance  with  Kolmogorov’s  2/3 rd  law. 


3.2.3  Fitting  the  representation  to  DNS  correlations 


To  evaluate  the  accuracy  of  the  anisotropy  approximation  developed  here,  we  first  evaluate  the  free 
parameters  in  the  model  by  fitting  to  correlations  determined  from  the  DNS  and  then  compare  the 
model  to  the  DNS  data.  The  fitting  procedure  is  described  in  [6]. 

The  fully  developed  turbulent  channel  flow  at  ReT  =  940  (Table  3)  is  used  for  the  fitting  and  eval¬ 
uation.  Correlations  were  computed  from  DNS  results,  the  details  of  which  can  be  found  in  [12]. 
This  is  a  highly  inhomogeneous  flow  and  so  the  extent  to  which  the  homogeneous  representation 
is  valid  will  also  be  examined. 

The  only  direction  of  inhomogeneity  is  the  one  perpendicular  to  the  channel  wall.  The  x  axis  is  in 
the  stream  wise  direction,  and  the  y  axis  is  normal  to  the  wall.  The  correlation  can  be  written  as 
Rij(y,  r)  =  (-u'(x)w' ;(x  +  r))  \X2=y.  The  homogeneous  model  does  not  satisfy  the  inhomogeneous 
continuity  equationor  the  symmetry  of  the  inhomogeneous  correlationexactly.  Instead,  we  will 
assume  that  Rij(y,  r)  is  locally  homogeneous  at  y.  This  of  course  implies  that  all  our  fitting 
parameters  are  functions  of  y. 

It  is  known  that  the  inertial  range  of  Rij(y,  r)  in  the  log  layer  is  self-similar,  with  a  similarity 
variable  r/y  [12,  28].  To  acknowledge  this  fact,  the  fitting  volume  V/  =  (r;  re(y)  <  r  <  rL(y)} 
is  defined  such  that  r^y)  is  proportional  to  y  near  the  wall  (Fig.  9).  Fig.  10  shows  that  the  relative 
error  in  the  model  correlation  is  less  than  20%  throughout  the  channel.  Also  (Fig.  10),  the  isotropic 
component  R/(r)  contributes  the  most,  followed  by  R/'(r  ),  Rd(r),  Rr- (r),  with  the  contribution 
from  R^(r)  being  negligible.  The  Kolmogorov  constant  6),.  calculated  from  (i\  has  a  value  quite 
close  to  2.0.  In  fact,  this  is  true  for  the  whole  channel  (Fig.  11),  where  Ck(y)  varies  between  1 .9  and 
2.2.  This  can  be  explained  from  Fig.  12,  where  {pb,  Pd ,  Pq}  can  be  seen  to  be  larger  than  p/  =  2/3 
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Figure  9:  Upper  and  lower  limits  of  the  range  of  r  over  which  error  minimization  is  performed, 
for  different  y/h  locations.  The  range  is  given  by  re  <  r  <  rL,  where  re(y)  =  (5C/;.)3'/4//(y)  and 
rL(y)  =  min(y/2,  0.1) 


throughout  the  channel,  satisfying  the  assumption  of  small-scale  isotropy.  Also  consistent  with 
this  assumption  is  the  fact  that  et]  is  nearly  isotropic. 

The  power-law  exponents  pb(y)  and  pdilj)  are  nearly  constant  throughout  the  channel  and  could 
be  approximated  as  Pb(y)  ~  1.4  and  pd(y )  ~  1  (Fig.  12).  The  power-law  exponent  pq(p)  can 
be  approximated  as  pg(y)  «  2.0  up  to  y/h  =  0.5  and  then  it  starts  to  fluctuate  near  the  center 
of  the  channel.  This  occurs  because  the  contribution  from  is  so  small  near  the  center  of  the 
channel  that  the  total  error  ||RMOD  —  RDNS||  is  insensitive  to  large  variations  in  j>q  .  Qualitatively, 
(Fig.  13)  the  model  fits  the  data  quite  well  for  the  normal  components,  when  the  reference  location 
(: y+  =  114)  is  in  the  log  layer  (Fig.  13).  The  inclination  of  principle  axes  are  well  represented  and 
the  magnitude  of  the  contour  levels  match.  The  main  shortcoming  of  the  model  lies  in  its  inability 
to  capture  the  effect  of  inhomogeneity,  which  is  most  obvious  in  the  component,  where  the 
isocontours  of  the  DNS  extend  out  further  in  the  positive  r2  half  of  the  plane  as  compared  to  the 
negative  r2  half,  due  to  the  presence  of  the  wall.  The  model  on  the  other  hand  is  symmetric  in  r,  and 
is  not  able  to  capture  this  aspect  of  the  data.  This  points  to  the  need  for  modeling  inhomogeneity 
in  the  representation. 


3.2.4  Anisotropy  due  to  Inhomogeneity 


The  results  of  fitting  the  homogeneous  model  to  the  inhomogeneous  correlations  from  the  channel 
flow  (figure  13),  clearly  show  the  shortcomings  on  the  homogeneous  model  in  representing  the 
cross-correlation  RV2.  Including  inhomogeneity  introduces  significant  complexity  in  the  represen¬ 
tation.  One  approach  to  this  more  complex  formulation  is  described  briefly  below. 

The  inhomogenous  two  point  correlation  Rij(x.,  x')  =  depends  on  two  position 

vectors  x  and  x',  and  can  therefore  also  be  mode  to  depend  on  any  two  linear  combinations  of 
these  two  position  vectors.  For  instance,  we  could  define  RT(y,  r)  =  (uj(x)uj(x')),  where  y  = 
r/x  +  (1  —  r/)x'  and  r  =  x'  —  x.  RtJ  satisfies  the  symmetry  /?7J(x.  x')  =  Rji(x.',  x).  It  turns 
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Figure  10:  Fractional  error  of  the  model  correlation  with  respect  to  the  DNS  correlation  at  various 
y/h  locations  across  the  channel.  The  bold  curve  shows  the  error  of  the  total  model  correlation. 
The  rest  of  the  curves  show  the  errors  of  the  truncated  model  correlation  at  different  y/ h  locations, 
with  various  components  subtracted  from  the  model  correlation  to  form  the  truncated  correlations. 


Figure  11:  Ck  obtained  from  cii(y) 
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Figure  12:  Power-law  exponents  {pb,Pd,PQ}  plotted  for  different  y/h  locations  in  the  channel.  pr 
has  been  fixed  at  2/3  for  the  whole  channel 


out  that  the  symmetry  has  the  simplest  form  for  rj  =  1/2,  i.e.  for  y  =  (x  +  x')/2,  we  have 
R}j2( y,  r)  =  R]{2( y,  —  r).  This  form  for  the  symmetry  is  easily  satisfied  by  (for  instance)  a  tensor 
that  depends  on  y,  is  even  in  r  and  is  symmetric  in  ij. 

The  representation  of  (ujfxjUj  (x,j)  is  therefore  best  done  in  the  mid-point  coordinates  (r/  =  1/2), 
because  the  symmetry  conditions  are  easier  to  satisfy.  From  here,  the  1/2  superscript  to  denote  the 
midpoint  coordinate  will  be  dropped. 


Two-point  correlation  in  the  mid-point  coordinate  is  Rij(y,  r)  =  (udy  —  r/2)uj(y  +  r/2)),  and 
the  two  continuity  conditions  that  need  to  be  satisfied  are: 
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A  form  that  satisties  continuity  can  then  be  given  as  follows: 


(46) 

(47) 


Rij{ y,  r)  =  eukejm^  <9+  [Hkn( r,  t)a(y)]  (48) 

where  tl3  is  a  tensor  that  satisfies  tkk  =  0.  The  symmetry  condition  can  be  easily  satisfied  by 
assuming  Hij( r,  t)  =  Hji(— r,  t).  We  also  implicitly  assume  that  does  not  vary  with  y.  That  is, 
the  relative  strengths  of  the  tensor  components  do  not  change  with  location.  This  is  in  fact  a  valid 
assumption  for  the  anisotropies  of  the  componentality  and  dimensionality  tensors  bij  and  in  the 
log  layer,  where  the  correlation  is  self-similar.  In  this  model,  the  inhomogenieties  of  Bl3  and  Dtj 
are  then  being  manifested  through  the  y  dependence  of  Btl(y)  and  Dtt  (  y  )  (which  are  not  the  same 
for  inhomogenous  turbulence). 
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Figure  13:  Isocontours  of  Rf^s(y,  r)  (DNS)  and  Rij(y ,  r)  (Model)  at  y+  =  114,  plotted  for  rz  =  0. 
The  correlations  in  the  left  column  were  calculated  from  DNS  of  channel  flow,  and  the  correlation 
in  the  right  column  from  the  best  fit  of  the  model  with  DNS.  The  contour  levels  for  a  given  corre¬ 
lation  component  (i.e.  for  the  same  row)  have  the  same  range. 
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Given  tij  =  tji,  four  invariant  forms  for  H%3  can  be  used: 

Hlj{  r,t)  =  fi(r)tij  (49) 

t)  =  f2{r)Sij{T  ■  t  ■  r)  (50) 

Hfj(  r,t)  =  /aWnr^r-t-r)  (51) 

=  /4(r)[ri(r  ■  t)j  +  r,-(r  ■  t)i]  (52) 

Without  any  loss  in  generality,  we  also  assume  that  the  origin  of  y  is  at  (x  +  x')/2,  and  model 
a°  (y  )  (which  multiplies  the  corresponding  //'))  as  a  quadratic  function  in  y,  i.e. 

aa(y)  —  ca  +  \%yk  +  Aastysyt  (53) 

Note  that  since  we  are  including  all  inhomogenous  effects  in  a“(y),  if  ca  =  /3c7,  then  it  implies 
that  =  /3X l  and  A°t  =  (3Ajt.  That  is,  any  constraints  amongst  aQ(y)  will  have  to  hold  across 
its  Taylor  series  expansion.  Therefore,  the  inhomogenous  representation  of  Rl3  for  a  given  t,3  is 
given  by 


Rij( y,r)  =J2eiikejmndl  d+  [H%n(r,  t)aQ(y)] 


(54) 


Expanding  d+  and  d  : 
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Rij( 0,  r)  can  therefore  be  written  as  a  sum  of  3  terms: 

Rij(®i r)  —  +  yW  +  2|(r) 
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where 
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Here,  R^fa  are  homogenous  tensors,  symmetric  in  r  and  ij,  Rf-a  arises  from  the  inhomogeniety 
(it  has  A"  in  every  term,  arising  out  of  first  derivative  of  a"(y)),  is  antisymmetric  in  r  and  ij  and 
R,y  also  arises  from  the  inhomogeniety  (due  to  Af)  and  is  symmetric  in  r  and  ij. 

If  all  Rfj'a( r)  have  the  same  power  law  (i.e.  fi  =  rp,  f2  =  rp~2,  f3  =  rp~A  and  /4  =  rp~  2), 

then  it  can  be  shown  that  {R'7-"}  are  linearly  dependent,  and  spanned  by  2  tensors.  However,  it 
can  also  be  shown  that  all  four  tensors  in  {Rs’“}  will  remain  linearly  independent.  Therefore, 
while  satisfying  the  self-consistency  condition,  we  can  independently  apply  constraints  for  Rffa 
and  R.y,  even  though  they  are  both  even  functions  in  r. 

This  representation  was  used  to  fit  channel  flow  from  the  DNS  at  Rer  =  940  and  the  results  for 
the  normal  components  of  the  correlation  tensor  are  the  same.  However,  the  antisymmetric  part 
breaks  the  symmetry  of  the  RV2  component  as  shown  in  figure  14.  This  is  an  improvement  over 
the  results  shown  in  figure  13,  especially  for  small  separations. 


3.3  Generalization  and  Testing  of  Finite-Volume  Wall-Bounded  OLES 

The  application  of  OLES  to  general  wall-bounded  flows  requires  a  synthesis  of  the  wall-bounded 
OLES  modeling  pursued  previously  using  spectral  representations  [39,  5]  and  the  finite-volume 
formulation  developed  and  evaluated  for  non-wall-bounded  flows  (see  section  3.1  and  [42,  26]). 
There  were  a  number  of  unanticipated  complications  in  this  generalization  of  optimal  LES,  the  two 
most  important  of  which  will  be  discussed  here.  These  have  been  investigated  in  the  context  of 
statistical  data  obtained  from  direct  numerical  simulation  of  channel  flow  at  Rer  =  940  [12]. 


27 


R\2  (DNS)  R12  (Model) 


Figure  14:  Isocontours  of  Rf^s(y,  r)  (DNS)  and  Rij(y,  r)  (Model)  at  y+  =  114,  plotted  for  rz  =  0. 
The  correlations  on  the  left  were  calculated  from  DNS  of  channel  flow,  and  the  correlation  in  the 
right  column  from  the  best  fit  of  the  model  with  DNS.  The  contour  levels  have  the  same  range. 

3.3.1  The  Spanwise  Flux  of  Mean  Momentum 


As  with  the  finite  volume  formulation  described  in  section  3.1,  in  wall-bounded  flows,  the  con¬ 
vective  flux  of  momentum  through  the  faces  of  the  finite  volumes  is  estimated  using  the  sum 
of  a  quandratic  and  linear  terms,  where  the  linear  term  is  expected  to  represent  the  transfer  of 
energy  to  the  unresolved  scales.  Using  DNS  statistical  data,  anisotropic  OLES  optimal  models 
were  developed  on  minimally  sized  stencils  (1x1x2)  for  simplicty.  In  this  inhomogeneous  flow, 
the  quadratic  terms  include  the  production  of  turbulent  kinetic  energy  through  interaction  with  the 
mean,  in  addition  to  the  nonlinear  turbulent  cascade  process.  As  a  result,  the  quadratic  term  is  a 
net  producer  of  kinetic  energy  through  much  of  the  simulation  domain.  However,  it  was  observed 
that  the  quadratic  term  slightly  underestimated  the  rate  of  production  in  the  central  region  of  the 
channel.  The  model  is  guaranteed  to  correctly  represent  the  energy  transfers  ( a  priori)  so  the  linear 
term  must  make  up  for  the  missing  production  by  being  slightly  anti-dissipative  near  the  channel 
centerline.  This  is  demonstrated  in  figure  15.  The  result  when  this  model  is  used  in  an  LES  is  an 
unstable  simulation. 

Further  investigation  of  this  phenomenon  identified  an  unlikely  culprit,  the  spanwise  flux  of  mean 
streamwise  momentum  Uw'  ( U  is  the  mean  velocity  and  w'  is  the  fluctuating  spanwise  veloc¬ 
ity)  that  appears  in  the  equation  for  streamwise  velocity  fluctuations  u'.  This  is  counter-intuitive 
because  this  term  cancels  out  of  the  equation  for  (■ u /2).  However,  the  OLES  models  ensure  correct 
a  priori  values  of  many  statistical  quantities,  including  those  involving  this  spanwise  flux.  In  par¬ 
ticular,  the  contribution  of  the  quadratic  term  to  the  spanwise  flux  U w'  was  underpredicted  by  the 
model  and  the  linear  term  had  to  make  up  the  difference.  However,  the  form  of  the  linear  term, 
which  involved  only  streamwise  velocities  was  incompatible  with  the  actual  term  U w',  which  is 
linear  in  the  (fluctuating)  spanwise  velocity.  To  correct  this  shortcoming,  the  model  was  general¬ 
ized  to  include  a  linear  in  w'  term  in  them  streamise  momentum  equations,  and  this  resulted  in  a 
stable  simulation.  The  near- wall  velocity  variances  from  this  simulation  are  shown  in  figure  16. 
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Figure  15:  Streamwise  energy  production  and  destruction  through  the  uw  term  by  the  LES  con¬ 
vective  terms  in  the  channel  as  determined  from  the  filtered  DNS,  the  quadratic  model  term  and 
the  linear  model  term. 


Figure  16:  Velocity  variances  from  the  filtered  DNS  (fDNS)  and  the  optimal  LES  (LES). 
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Figure  17:  Contribution  of  the  convection,  viscous  and  pressure  terms  to  the  evolution  of  the 
filtered  velocity  variances  as  determined  from  the  filtered  DNS  and  optimal  LES  model. 


3.3.2  Near-wall  Pressure  Contributions 


The  obvious  shortcoming  of  the  LES  results  shown  in  figure  16  is  that  the  streamwise  fluctuations 
are  over-predicted  compared  to  the  filtered  DNS,  and  the  other  velocity  fluctuation  components  are 
somewhat  underpredicted.  A  likely  source  of  this  discrepency  is  shown  in  figure  17  in  which  the 
agregated  terms  in  the  evolution  equations  for  ( u /2)  and  (v'2)  are  shown.  Notice  that  the  primary 
discrepency  between  the  LES  model  terms  and  the  DNS  terms  is  in  the  pressure  contribution,  which 
is  responsible  for  the  transfer  of  energy  from  the  streamwise  to  the  other  velocity  components. 

In  this  OLES  model  formulation,  the  pressure  is  determined  by  imposing  a  discrete  divergence 
free-constraint  on  the  volume  averaged  velocities,  as  in  a  standard  staggered  grid  finite  volume 
scheme.  This  is  a  modeling  ansatz  is  based  on  the  observations  of  [20]  that  in  finite  volume  LES 
of  isotropic  turbulence  this  is  very  nearly  the  the  optimal  model.  The  presence  of  the  wall  appears 
to  be  upsetting  that  conclusion.  There  are  several  ways  to  interpret  the  discrepency  in  the  pressure 
term.  The  most  useful  is  to  observe  that  the  magnitude  of  the  discrete  divergence  of  the  convective 
term  is  greatly  underestimated  by  the  model  relative  to  the  filtered  DNS  (figure  18),  resulting  in  a 
lower  pressure  contribution.  This  is  a  property  of  the  convective  terms  that  the  OLES  formulation 
does  not  control  for.  The  proposed  solution  to  this  problem  is  to  reformulate  the  OLES  model  to 
estimate  the  (discrete)  divergence-free  projection  of  the  convection  term.  This  proposal  is  currently 
being  evaluated. 


3.4  Conclusions 

3.4.1  Theory-based  Finite- Volume  OLES 


The  correlations  required  for  the  finite-volume  OLES  formulation  can  be  obtained  from  Kol¬ 
mogorov  inertial-range  theory,  small-scale  isotropy,  and  the  quasi-normal  approximation.  Further, 
LES  models  resulting  from  these  correlations  perform  well.  These  approximations  will  be  valid 
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Figure  18:  Variance  of  the  divergence  of  the  LES  convection  term  N,  as  determined  from  the  DNS 
(iV)and  the  optimal  LES  model  (Nm). 


provided  the  turbulence  at  the  filter  scale  and  smaller  is  locally  homogeneous  and  isotropic.  Except 
for  very  near  walls  and  other  strong  inhomogeneities,  this  is  expected  to  be  a  good  approximation 
provided  the  Reynolds  number  is  sufficiently  high  and  the  filter  width  is  sufficiently  small. 

The  finite-volume  OLES  models  developed  here  are  expressed  as  discrete  quadratic  and  linear 
operators  that  represent  the  convective  momentum  flux.  These  operators  depend  explicitly  on 
two  flow  dependent  parameters:  the  dissipation  rate  e  and  the  nondimensional  parameter  7  = 
A e/u3  (the  ratio  of  the  grid  size  A  to  the  large  turbulence  scale).  The  modeling  approach  allows 
the  parameters  e  and  u2  to  be  determined  dynamically  and  accurately  in  a  running  LES.  But  it 
was  also  found  that  asymptotic  operators  for  7  — >  0  yield  results  consistent  with  the  finite  7 
operators  with  7  as  large  as  0.08.  In  most  cases,  therefore,  the  asymptotic  operators  are  sufficiently 
accurate,  meaning  that  only  e  needs  be  determined  dynamically.  The  process  for  determining  e  is 
considerably  simpler  than  the  usual  dynamic  procedure  [13],  as  no  test  filter  is  required. 

The  quadratic  and  linear  operators  arising  from  the  LES  optimization  are  broadly  similar  to  com¬ 
mon  finite- volume  operators.  However,  there  are  significant  quantitative  differences  between  the 
OLES  operators  and  standard  finite-volume  operators,  which  should  be  interpreted  as  part  of  the 
model  for  the  effects  of  subgrid  scales.  The  linear  part  of  the  operators,  which  play  the  role  of 
the  subgrid  stress  model  by  dissipating  resolved  energy,  are  different  for  the  normal  and  tangential 
velocity  components  in  any  grid  direction,  with  the  tangential  component  dissipation  significantly 
lower  at  high  wavenumbers  than  a  standard  fourth-order  approximation  to  the  second  derivative. 
This  is  one  of  the  features  of  the  OLES  approach;  the  spectral  distribution  of  subgrid  dissipation 
is  tailored  to  be  consistent  with  the  statistics  of  turbulence.  Finally,  note  that  the  quadratic  OLES 
operators  on  average  produce  resolved- scale  energy  and  the  linear  operators  dissipate  energy.  The 
dissipation  of  the  linear  term  is  determined  to  ensure  that  the  combined  contributions  of  quadratic 
and  linear  terms  to  the  resolved  energy  evolution  match  the  required  total  dissipation. 

These  results  indicate  that  practical  finite-volume  OLES  models  can  be  formulated  without  DNS 
correlation  inputs,  at  least  for  the  circumstances  for  which  most  LES  models  are  designed  (small- 
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scale  isotropy).  The  advantage  of  the  current  formulation  is  that  it  allows  the  effects  of  subgrid 
turbulence  and  numerical  discretization  to  be  treated  in  a  unified  way  and  consistent  with  the  statis¬ 
tical  properties  of  turbulence.  It  is  both  remarkable  and  encouraging  that  the  simple  considerations 
employed  here  are  sufficient  for  this  modeling. 

Several  immediate  generalizations  are  possible  using  the  isotropic  statistical  models  introduced 
here.  These  include  LES  of  scalar  transport  and  modeling  on  anisotropic  and  inhomogeneous  grids, 
as  one  would  encounter  in  most  applications.  Also  important  is  the  generalization  to  anisotropic 
correlations  for  LES  models  applicable  when  isotropy  at  the  grid  scale  is  not  a  valid  assumption. 
This  occurs,  for  example,  in  near-wall  turbulence.  Isotropy  places  a  strong  constraint  on  the  form 
of  the  multi-point  correlations.  Without  isotropy,  even  the  two-point  second-order  correlation  is 
difficult  to  represent  (see  [6]  and  the  references  therein).  The  possibility  raised  in  section  3.1.6 
that  the  I3  contribution  to  the  OLES  estimation  equations  might  be  neglected  or  crudely  modeled 
is  thus  particularly  important,  since  the  underlying  three-point  third-order  correlation  is  the  most 
difficult  to  model.  However,  the  robustness  of  this  result  to  anisotropy  and  inhomogeneity  of  the 
turbulence  and  the  LES  grid  is  far  from  clear. 


3.4.2  Modeling  Anisotropy  of  Two-point  Correlations 


The  results  show  that  the  homogeneous  model  fits  the  correlations  calculated  from  DNS  data  rather 
well,  both  quantitatively  and  qualitatively,  though  there  are  shortcomings  to  the  fit,  arising  from 
limitations  of  the  model.  In  particular,  to  improve  the  representation  of  the  small  separation  cor¬ 
relations  discussed  here  for  flows  that  are  inhomogeneous,  the  effects  of  the  inhomogeneity  on 
the  correlation  needs  to  be  represented,  particularly  the  breaking  of  the  symmetry  expressed  as 
Rij( r)  =  Rij(— r).  A  preliminary  inhomogeneous  formulation  was  proposed  and  tested,  with 
encouraging  results. 

The  primary  advantage  of  models  of  the  type  developed  here  for  the  small  separation  correlation  is 
that  they  have  a  finite  number  of  parameters  (17  including  the  power-law  exponents  for  the  homo¬ 
geneous  model).  If  we  neglect  the  part  of  the  correlation  depending  on  Q*jk  (which  accounts  for 
<  5%  of  the  correlation),  then  our  model  can  be  expressed  in  terms  of  just  9  free  parameters.  A 
similar  SO(3)  (up  to  /  =  2)  decomposition  for  Rrj (r)  that  has  even  parity  in  r,  is  symmetric  in  ij 
and  satisfies  the  homogenous  continuity  equation  will  have  8  free  parameters  [2]  (after  fixing  the 
isotropic  power  law)  -  thus,  the  number  of  DOFs  in  our  model  is  comparable  to  that  of  the  simplest 
possible  anisotropic  SO(3)  representation.  At  any  rate,  a  finite  number  of  measurements  of  the  cor¬ 
relation  will  be  sufficient  to  parametrize  our  model,  effectively  reconstructing  the  correlation  from 
a  small  amount  of  data.  This  is  useful  in  interpreting  experimental  data,  but  perhaps  more  impor¬ 
tant,  it  is  useful  for  reconstructing  the  turbulence  correlations  from  correlations  computed  from  a 
large  eddy  simulation.  Applying  the  LES  filter  to  the  correlation  model  yields  a  representation  of 
the  correlation  of  the  LES  velocities,  which  can  be  fit  to  the  LES  data  similar  to  the  fits  to  DNS 
performed  in  this  paper.  The  parameters  so  determined,  when  applied  to  the  unfiltered  correlation 
representation,  define  the  turbulence  correlation  consistent  with  the  correlations  of  the  LES  fields. 

Seventeen  is  a  large  number  of  parameters,  which  raises  the  questions  as  to  how  practical  it  might 
be  to  parametrize  such  a  model  using  LES  (or  experimental)  data.  However,  this  is  a  second 
rank  correlation  tensor.  Correlations  for  just  two  non-zero  values  of  r  yield  18  data  items.  If  the 
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correlation  (or  the  filtered  correlation)  is  sampled  on  a  rectangular  grid  (e.g.  an  LES  grid),  with 
nearest  neighbor  separations  (7  points  total)  that  would  amount  to  33  data  items  (accounting  for 
symmetries),  and  sampling  on  a  cubical  grid  of  27  separations  yields  123  independent  data  items. 
Clearly,  data  from  a  reasonable  number  of  small  separations  should  be  sufficient  to  parametrize  the 
model.  An  effort  to  extract  Rij  in  this  manner  from  LES  statistics  was  undertaken  by  Bhattacharya 
[4],  where  encouraging  results  were  obtained  for  extracting  the  normal  Reynolds  stress  components 
Baa  from  LES  correlations. 

While  the  results  described  above  are  encouraging  for  the  representation  of  anisotropy  in  the  two- 
point  second-order  velocity  correlation,  the  optimal  LES  formulation  also  needs  a  model  for  the 
three-point  third-order  correlation.  The  complexity  of  such  a  model  will  be  much  greater  than  for 
the  two-point  correlation  discussed  here,  so  much  so  that  this  approach  to  anisotropy  representation 
of  the  correlations  underlying  the  LES  formulation  now  appear  impractical.  The  current  results  for 
the  two-point  correlation  will  be  useful  for  reconstructing  second  order  statistics  from  LES  fields, 
but  another  approach  to  representing  the  anisotropy  in  optimal  LES  models  appears  to  be  needed. 

A  much  more  direct  and  simplier  approach  is  currently  being  pursued.  It  is  based  on  the  modeling 
ansatz  that  the  anisotropy  of  the  interaction  between  rsolved  and  unresolved  scales  can  be  char¬ 
acterized  by  anisotropy  in  the  dissipation  inferred  from  an  LES  using  the  techniques  described  in 
[26].  This  approach  is  currently  being  evaluated  for  its  effectiveness. 


3.4.3  Finite  Volume  Wall-bounded  OLES 


The  execution  of  finite  volume  OLES  for  wall-bounded  turbulence  encountered  several  unantic¬ 
ipated  roadblocks,  which  have  been  or  are  being  overcome.  As  with  previous  development  of 
OLES  models,  the  resolution  of  the  problems  encountered  was  facilitated  by  the  formal  structure 
of  the  modeling  approach,  which  provides  a  theoretical  framework  in  which  the  expectations  for 
the  statistical  properties  of  the  model  are  analyzed.  The  analysis  also  provides  insight  with  poten¬ 
tial  application  for  wall-bounded  LES  in  general.  In  particular,  one  of  the  common  shortcomings 
of  wall-bounded  LES,  especially  with  unresolved  wall  layers,  is  the  over  prediction  of  streamwise 
velocity  fluctuations,  just  as  observed  here.  It  may  be  that  this  shortcoming  is  commonly  related 
to  inadequacies  in  the  treatment  of  pressure  and  the  continuity  constraint,  as  is  the  case  here. 
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